do_RSA_v12

Load libraries and define parameters

NB: this Rmd can be used for all types of ratings: emotion, arousal, valence

Abbreviations: - sub : subjects(s) - rat : rating(s)

Load libraries

library(tidyverse)
library(future)
library(furrr)
library(tictoc)
library(RNifti)
library(proxy) # distances
library(profvis)
library(DT)
library(formattable)
library(janitor)
library(ComplexHeatmap)
library(circlize)  # For color mapping functions
library(grid)  # For gpar()
# library(pheatmap)
# library(heatmaply)
library(papayaWidget)
library(ggstatsplot)


source("funs_V12_ROIs.R")

Define paths and rsa_[flavour] specific option

bd="/data00/leonardo/RSA/analyses" # change to your path
bd_ratings = paste0(bd,"/RATINGS")
bd_atlases = paste0(bd,"/ROIS_REPO")

# Directory to store result for rsa explorer app. Must be created manually
bd_results = paste0(bd,"/RSA_ROI_APP/results_RSA_ROI") 


ratings_type <- c("emotion","arousal","valence","aroval")

rsa_flavour="rsa_ROI"

# Do you want to remove neutral cope(s) from the RSA analysis?
remove_neutrals = "YES"

# choose the type of betas to be used, e.g.
# - one_ev_per_movie [_minus_neutral]
# - emotion_predictors [_minus_neutral]
copes_type <- "one_ev_per_movie"


# Choose one atlas, e.g.
# Schaefer100.nii.gz
# Schaefer200.nii.gz
# Yeo7.nii.gz
atlas_filename <- "Yeo7.nii.gz"
# Atlases repos:
# https://www.lead-dbs.org/helpsupport/knowledge-base/atlasesresources/cortical-atlas-parcellations-mni-space/
# https://github.com/neurodata/neuroparc
# https://www.fmrib.ox.ac.uk/datasets/brainmap+rsns/

# Choose the distance metric to use for fmri RDM
# The metric for ratings RDM *should be* euclidean, since arousal and valence
# ratings have only one value
# Supported methods: 'pearson', 'spearman', 'euclidean', 'cosine' or 'mahalanobis'
dist_method_rating <- "euclidean"
dist_method_fmri = "euclidean"
dist_method_rsa = "pearson"

# Vector of zeropadded sub_ids
subs_file <- "/data00/leonardo/RSA/sub_list.txt"
subs <- sprintf("%02d",readLines(subs_file) %>% as.numeric)

# ---------- ONLY TOP RATERS BELOW ------------
subs <- c("02","03","12","11","22","09","29","28","26","32","23","15","20","19")

Necessary columns for ratings/models

# df <- read_csv(paste0(bd_ratings,"/emotion_ratings.csv"))
# df %>%
#   select(sub, video, r_happy, r_fear, r_pain, movie_cope_number)
# 
# 
# # to remove specific *rows*
# videos_to_remove <- c("AK_A_32.mp4", "AK_A_04.mp4")
# 
# df %>%
#   filter(sub == "02") %>% 
#   filter(!video %in% videos_to_remove)  

Create copes metadata (to identify and remove neutrals)

Before anything else, we need to identify the neutral copes and their associated ratings.

To this aim we load the emotion ratings which contain also the correspondence between movie, cope and rating. Then we select only one sub (the order is the same for all subs)

metadata_copes <- read_csv(paste0(bd_ratings, "/emotion_ratings.csv")) %>% 
  filter(sub == subs[1]) %>% 
  select(high_low_code, movie_cope_number, emotion) %>% 
  rename(cope = movie_cope_number, label = high_low_code)
  
metadata_copes

## # A tibble: 56 × 3
##    label      cope emotion
##    <chr>     <dbl> <chr>  
##  1 AK_A_high     1 anger  
##  2 AK_A_low      2 anger  
##  3 AK_D_high     3 disgust
##  4 AK_D_low      4 disgust
##  5 AK_F_high     5 fear   
##  6 AK_F_low      6 fear   
##  7 AK_H_high     7 happy  
##  8 AK_H_low      8 happy  
##  9 AK_N_high     9 neutral
## 10 AK_N_low     10 neutral
## # ℹ 46 more rows

Create a df_path_copes with the location of the 56 cope niis from the one_movie_per_ev model

Extract the pathname of all copes using the list.files() function. Also define a copes_numba vector with all the copes numbers.

NB: The cope numbers in the cope column are NOT zeropadded since this is how they come out from FSL Feat

# load df_path_copes and add metadata (label and emotion/neutral)
df_path_copes <- import_df_path_copes(bd, copes_type, rsa_flavour) %>% 
  inner_join(metadata_copes, by="cope")

# If remove_neutrals is "YES", remove neutral copes and arrange by sub/cope (just to be sure) 
if (remove_neutrals == "YES") {
  df_path_copes <- df_path_copes %>%  
  filter(emotion != "neutral") %>% 
  arrange(sub,cope)
} 

copes_numba <- df_path_copes$cope %>% unique

cat(
  "there are",length(copes_numba)," copes including ",
  nrow(df_path_copes %>% filter(sub==subs[1], emotion=="neutral"))," neutral copes \n"
)

## there are 48  copes including  0  neutral copes

# df_path_copes
# length(copes_numba)

Read atlas and create a vector of atlas labels

atlas_path <- paste0(bd_atlases,"/",atlas_filename)

atlas_nii <- readNifti(paste0(bd_atlases,"/",atlas_filename))
region_labels <- atlas_nii[atlas_nii > 0] %>% unique %>% sort

Calculate RATINGS RDMs

# NB: ratings_type is now a *vector*, e.g.: 
# ratings_type <- c("emotion","arousal","valence","aroval")
# to test the fn below for one rating type : do_RDMs_rats("emotion")

do_RDMs_rats <- function(ratings_type, remove_neutrals) {

  ratings_path <- paste0(bd_ratings,"/",ratings_type,"_ratings.csv")
  rats <- read_csv(ratings_path)
  
  if (remove_neutrals == "YES") {
    rats <- rats %>% filter(emotion != "neutral")
  }
  
  rdm <- rats %>%
    filter(sub %in% subs) %>%
    select(sub, starts_with("r_")) %>%
    group_by(sub) %>%
    nest %>%
    mutate(!!paste0("RDM_",ratings_type) := data %>% map(~ DDOS_vec(.x, dist_method_rating))) %>%
    ungroup %>% 
    # remove the sub with the prospect of merging rdms of different ratings/models
    select(!data)
  
  return(rdm)
  
}


# calculate rdms for all ratings_type and return one column for each ratings_type
RDMs_rats <- map_dfc(ratings_type, ~ do_RDMs_rats(.x, remove_neutrals)) %>%
  # remove the duplicated sub_ columns (generated by do_RDMs_rats) 
  janitor::clean_names() %>% 
  mutate(sub = sub_1) %>% 
  select(sub, starts_with("rdm"))


# # example heatmap
# plot_tril(
#   RDMs_rats[1,]$rdm_emotion, model = "emotion",
#   fontsize=7, side_mm=100, reord = FALSE
# )

Calculate FMRI RDMs

IMPORTANT : Since we already removed neutral from df_path_copes we don’t need to filter anything anymore here.

# Fn to calculate the rdm_fmri for one sub
# Test for one sub: do_RDM_fmri(subs[1], copes_numba, df_path_copes)

do_RDM_fmri <- function(sub_id, copes_numba, df_path_copes) {
  
  # load all the copes nii.gz into a df
  df_copes <- load_sub_copes(sub_id, copes_numba, df_path_copes)
  
  # the code below does the calculation of the tril of RDM_fmri for all rois for one sub
  
  one_sub_RDM_fmri <- tibble(
    sub = sub_id,
    roi = region_labels
  ) %>%
  
    # extract idx_roi for atlas voxels in that roi
    mutate(idx_roi = roi %>% map(~ which(atlas_nii == .x)) ) %>%
    
    # extract the df_copes for the voxels in that region (idx_roi)
    mutate(df_copes_region = idx_roi %>% map(~ df_copes[.x,]) ) %>%
    
    # calculate tril of RDM_roi (output as numeric vector)
    mutate(rdm_fmri = df_copes_region %>% map(~ DDOS_vec(t(.x), dist_method_fmri)) ) %>% 
    
    # select only sub, roi (numba) and RDM_roi
    select(sub, roi, rdm_fmri) 
  
  return(one_sub_RDM_fmri)
}



# Calculate RDMs_fmri for all subs
plan(multisession, workers = 5)

RDMs_fmri <- subs %>% future_map_dfr(
  ~ {
    paste0("Calculating RDMs for sub ",.x,"\n") %>% cat
    do_RDM_fmri(.x, copes_numba, df_path_copes)
  }
)

## Calculating RDMs for sub 02
## Calculating RDMs for sub 03
## Calculating RDMs for sub 12
## Calculating RDMs for sub 11
## Calculating RDMs for sub 22
## Calculating RDMs for sub 09
## Calculating RDMs for sub 29
## Calculating RDMs for sub 28
## Calculating RDMs for sub 26
## Calculating RDMs for sub 32
## Calculating RDMs for sub 23
## Calculating RDMs for sub 15
## Calculating RDMs for sub 20
## Calculating RDMs for sub 19

plan(sequential)

# RDMs_fmri

Calculate RSA

# The dist_method_rsa is set above.
# can be pearson, spearman
# dist_method_rsa <- "spearman"

# Join fmri and rats RDMs and create a copy of the rat RDM for each sub/roi
RSA <- right_join(RDMs_fmri, RDMs_rats, by = "sub") %>%
  group_by(sub) %>%
  
  # Dynamically create rsa_ columns for every ratings type
  mutate(across(
    all_of(paste0("rdm_", ratings_type)),
    ~ map2_dbl(.x, rdm_fmri, ~ cor(.x, .y, method = dist_method_rsa)),
    .names = "rsa_{.col}"
  )) %>%
  mutate(across(starts_with("rsa_"), round, 3)) %>%
  ungroup()

## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `across(starts_with("rsa_"), round, 3)`.
## ℹ In group 1: `sub = "02"`.
## Caused by warning:
## ! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
## Supply arguments directly to `.fns` through an anonymous function instead.
## 
##   # Previously
##   across(a:b, mean, na.rm = TRUE)
## 
##   # Now
##   across(a:b, \(x) mean(x, na.rm = TRUE))

# Get the mean RSA for each ratings_type
RSA_mean <- RSA %>%
  select(roi, starts_with("rsa_")) %>%
  group_by(roi) %>%
  reframe(
    across(
      starts_with("rsa_"),
      list(mean = ~ round(mean(.x, na.rm = TRUE), 2))
    )
  ) %>% 
  rename_with(
    ~ str_replace(.x, "rsa_rdm_", ""),
    starts_with("rsa_")
  )

RSA_mean %>% datatable()

Calculate mean tril for ratings and for fmri (for each roi in the fmri case)

# Takes a df and a "column" containing a nested vector in each row
# and returns the element-wise mean of that vector across rows
get_elementwise_mean <- function(mydf, nested_col, fun=mean) {
  m <- matrix(
    mydf[[nested_col]] %>% unlist,
    nrow = nrow(mydf),
    ncol = length(mydf[[nested_col]][[1]]),
    byrow = TRUE
  )
  
  m_summary <- apply(m, 2, fun)
  return(m_summary)
}


# mean RDMs_rats
mean_RDMs_rats <- RDMs_rats %>% select(starts_with("rdm_")) %>% colnames %>%  
  map_dfc(~{
    mean_tril <- get_elementwise_mean(RDMs_rats, .x, fun = mean)
    tibble(!!.x := mean_tril)
  })



# mean RDMs_fmri
plan(multisession, workers = 5)

mean_RDMs_fmri <- RDMs_fmri$roi %>% 
  future_map_dfr(~{
    mean_tril <- get_elementwise_mean(RDMs_fmri %>% filter(roi == .x), "rdm_fmri")
    tibble(roi = .x, mean_tril = mean_tril)
  }) %>% 
  group_by(roi) %>% 
  nest

plan(sequential)

Plot mean RDMs_rats

# Plot mean RDMs_rats
ratings_type %>% 
  map(
    ~ plot_tril(
      mean_RDMs_rats[[paste0("rdm_",.x)]], model = .x,
      fontsize=8, side_mm=200, reord = TRUE
    )
  )

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

Plot mean RDMs_fmri for all ROIs

mean_RDMs_fmri

## # A tibble: 7 × 2
## # Groups:   roi [7]
##     roi data                 
##   <dbl> <list>               
## 1     1 <tibble [15,792 × 1]>
## 2     2 <tibble [15,792 × 1]>
## 3     3 <tibble [15,792 × 1]>
## 4     4 <tibble [15,792 × 1]>
## 5     5 <tibble [15,792 × 1]>
## 6     6 <tibble [15,792 × 1]>
## 7     7 <tibble [15,792 × 1]>

seq(nrow(mean_RDMs_fmri)) %>% 
  map(
    ~ plot_tril(
      mean_RDMs_fmri[.x,]$data, model = paste0("fmri_roi_",.x),
      fontsize=8, side_mm=200, reord = FALSE
    )
  )

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

Elements for the app

RSA table

RSA_mean %>% datatable()

Boxplot for one roi

roi_numba = 6

RSA %>% 
  select(sub, roi, starts_with("rsa_")) %>% 
  filter(roi == roi_numba) %>% 
  select(starts_with("rsa_")) %>% 
  pivot_longer(cols = starts_with("rsa_"), names_to = "model") %>%
  mutate(model = str_replace(model, "rsa_rdm_","")) %>% 
  ggwithinstats(
    x = model, y = value,
    type = "p",   # p, np, r, b
    results.subtitle = TRUE
  )

View one ROI

prepare_papaya <- function(atlas_filename, roi_numba) {
  
  f <- function(nii_filename) {
    return(paste0(bd_atlases,"/", nii_filename))
  }
  
  nii_atlas <- RNifti::readNifti( f(atlas_filename) )
  # view(nii_atlas)
  
  nii_atlas_thr <- ifelse(nii_atlas == roi_numba, 1, 0)
  RNifti::writeNifti(nii_atlas_thr, template = nii_atlas, paste0(bd_atlases,"/tmp_ROI.nii.gz"))
  
  papaya(
    c(f("Dummy.nii.gz"), f("MNI152_T1_2mm_brain.nii.gz"), f("tmp_ROI.nii.gz")),
    options = list(
      papayaOptions(alpha = 1, lut = "Grayscale"),
      papayaOptions(alpha = 0.5, lut = "Grayscale"),
      papayaOptions(alpha = 0.5, lut = "Red Overlay", min = 0, max = 5)
    ),
    interpolation = FALSE,
    orthogonal = FALSE
  )
  
}

prepare_papaya(atlas_filename, roi_numba = 6)

Plot MEAN tril for fmri and all models for one roi

global_reord = FALSE
global_fontsize = 7
global_side_mm = 100

roi_numba = 6

# fmri one roi
tril_fmri <- plot_tril(
  mean_RDMs_fmri[roi_numba,]$data, model = paste0("fmri_roi",roi_numba),
  fontsize = global_fontsize, side_mm = global_side_mm, reord = global_reord
)


# emotion
tril_emotion <- plot_tril(
  mean_RDMs_rats$rdm_emotion, model = "emotion",
  fontsize = global_fontsize, side_mm = global_side_mm, reord = global_reord
)

# arousal
tril_arousal <- plot_tril(
  mean_RDMs_rats$rdm_arousal, model = "arousal",
  fontsize = global_fontsize, side_mm = global_side_mm, reord = global_reord
)

tril_fmri + tril_emotion

Save results in RSA_ROI_APP

Make sure you already created the directory /data00/leonardo/RSA/analyses/RSA_ROI_APP/results_RSA_ROI

# We save the following into /data00/leonardo/RSA/analyses/RSA_ROI_APP
# RSA
# RSA_mean
# atlas_filename
# mean_RDMs_fmri
# mean_RDMs_rats
# df_path_copes
# subs


create_filename_results <- function(atlas_filename, copes_type, remove_neutrals) {
  atlas_root <- str_replace(atlas_filename, ".nii.gz", "")

  if(str_detect(copes_type,"one_ev_per_movie")) {EV_abbr <- "_OEPM"}
  if(str_detect(copes_type,"emotion_predictors")) {EV_abbr <- "_EP"}
  if(str_detect(copes_type,"emotion_high_low_predictors")) {EV_abbr <- "_EHLP"}
  
  option_MN <- ifelse(str_detect(copes_type,"_minus_neutral"),"_MN","")

  option_NN <- ifelse(remove_neutrals == "YES", "_RN", "")
  
  results_filename <- paste0(atlas_root, EV_abbr, option_MN, option_NN)
  return(paste0(bd_results,"/",results_filename,".RData"))
  
}

results_file <- create_filename_results(atlas_filename, copes_type, remove_neutrals)
results_file

## [1] "/data00/leonardo/RSA/analyses/RSA_ROI_APP/results_RSA_ROI/Yeo7_OEPM_RN.RData"

# save the results
save(RSA, RSA_mean, atlas_filename, mean_RDMs_fmri, mean_RDMs_rats, df_path_copes, subs, file = results_file)

# load them again
# load(results_file)
# load("/data00/leonardo/RSA/analyses/RSA_ROI_APP/results_RSA_ROI/Yeo17_OEPM_MN_RN.RData")